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We propose an effective permittivity modei to homogenize an array of long tliin e-negative rods 
arranged in a periodic lattice, ft is proved that the effect of spatial dispersion in this electromagnetic 
crystal cannot be neglected, and that the medium supports dispersionless modes that guide the 
energy along the rod axes. It is suggested that this effect may be used to achieve sub-wavelength 
imaging at the infrared and optical domains. The reflection problem is studied in detail for the case 
in which the rods are parallel to the interfaces. Full wave numerical simulations demonstrate the 
validity and accuracy of the new model. 
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I. INTRODUCTION 

In recent years there has been a lot of interest in the 
propagation of electromagnetic waves in artificial mate- 
rials, and particularly in materials with negative index 
of refraction [1, 2]. The research has been mainly driven 
by the possibility of using these materials to miniatur- 
ize several devices and waveguides, and develop "per- 
fect lenses" able of focusing electromagnetic radiation 
with sub- wavelength resolution [3] . Recently, a different 
approach to achieve sub-wavelength resolution was pro- 
posed in [4] and demonstrated experimentally in [5]. In 
[5] the idea is to use an artificial material formed by an ar- 
ray of perfectly conducting wires (wire medium) to guide 
the waves from the input plane to the output plane, and 
then reconstruct the image " pixel by pixel" , exploring a 
Fabry-Perot resonance that in the wire medium occurs 
simultaneously for all the spatial harmonics. The reso- 
lution of this transmission device is only limited by the 
lattice constant, i.e. by the spacing between the wires. 
The problem with the configuration studied in [5] is that 
it is limited to the microwave regime because in the opti- 
cal domain perfect electric conducting materials are not 
available. Nevertheless, we will suggest in this paper that 
it may still be possible to use the same concept to achieve 
sub- wavelength resolution at the infrared and optical do- 
mains, provided e-negative (ENG) rods are used to guide 
the electromagnetic radiation instead of metallic wires. 
At the infrared and optical frequencies all metallic ma- 
terials have permittivity with negative real part. It is 
known that their dielectric constant is well represented 
by the Drude model e = 1 — iJ^juj^ , where lo^ is the 
plasma frequency (for simplicity the lossless model was 
considered). Thus, we envision that either silver or gold 
rods may used to fabricate transmission devices that are 
able to propagate the subwavelength information of an 
image. 

With this motivation, we will investigate in this paper 
the propagation of electromagnetic waves in an artificial 
medium formed by thin ENG rods (see Fig. 1), and we 
will show that the structure can be homogenized and de- 
scribed by an effective permittivity tensor provided spa- 



tial dispersion is taken into account. We will show that 
our homogenization model predicts that the rod medium 
may support a mode that propagates along the axes of 
the rods with the same phase velocity, independently of 
the transverse wave vector. As proved in [4], this is a 
key property to operate the material in the canalization 
regime and achieve sub-wavelength resolution. 




FIG. 1: Medium formed by long thin ENG rods arranged in 
a square lattice. 

Previous works related with the homogenization of 
ENG rods are discussed next. In [6] the oblique prop- 
agation of electromagnetic waves through an array of 
aligned fibres was examined, and the associated prob- 
lem of numerical instability discussed. In [7] a method 
was proposed to compute analytically the band struc- 
ture of wire mesh crystals. In [8] a model was proposed 
to homogenize a medium with metallic rods, but only the 
on-plane case was considered. In [9] an homogenization 
model was proposed to characterize a structure similar 
to the one studied in this paper. It was demonstrated 
that the artificial material was characterized by spatial 
dispersion, and that the band structure of the photonic 
crystal had several branches. However, the results of [9] 
are restricted to the case in which the permittivity of 
the rods follows a Drude-type plasma model, and besides 
that the derived formulas for the effective permittivity 
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arc very cumbersome, lead to non-analytical dispersion 
characteristics, and more importantly hide the physics 
of the problem. In this paper, we will derive a simpler 
and more intuitive model, describe new phenomena, and 
prove that the new approach can characterize accurately 
the electrodynamics of the artificial medium. 

The homogenization of the wire medium is also closely 
related with the subject under study. In fact, the case of 
perfectly conducting wires can be regarded as the limit 
situation in which the permittivity of the ENG rods is 
—00. In [10] an homogenization model was derived for 
the wire medium, and it was proved that this artificial 
material suffers from strong spatial dispersion even for 
very long wavelengths. Later, in [11, 12], these results 
were generalized for 2D- and 3D lattices of connected and 
unconnected wires. In [13, 14] the reflection problem in 
a wire medium slab was investigated, and in [15] it was 
proved that in general an additional boundary condition 
is necessary to determine the scattering parameters. 

The paper is organized as follows. Firstly, we will 
characterize the electric polarizability of a dielectric rod. 
Then, we will derive the homogenization rudiments nec- 
essary to calculate the effective permittivity of the struc- 
ture. In section IV, the electromagnetic modes supported 
by the rod medium are characterized. In section V, we 
study the reflection problem in a finite slab, assuming 
that the rods are parallel to the interfaces. The derived 
results are compared with full wave numerical simula- 
tions. Finally, in section VI the conclusions are pre- 
sented. 

In this work we assume that the fields are monochro- 
matic with time variation e~^^"*. 



II. POLARIZABILITY OF A DIELECTRIC ROD 

Let us consider a ENG rod with radius R and (rela- 
tive) permittivity e = e{ui). The rod is oriented along 
the z-direction. In this section, we calculate the com- 
ponent azz of the electric polarizability tensor (actually, 
since the rod is infinite along z, we are going to calculate 
the polarizability per unit of lenght). To this end, we 
consider that a plane wave with magnetic field along the 
y-direction illuminates the dielectric rod. The wave vec- 
tor of the incident wave is k = {k^, 0, kz), where k.k = 
and /? = w/c is the free-space wave number. The incident 
electric field along the z-direction is: 

£1^"'= = £;oe-^'=»'^e--'*^'^ (1) 

OC 

= Eoe-^'''^ ^ j|^l(fc^ ^r)e^"^ 

n= — oo 

where Jn is the J-Bessel Function of first kind and order 
n, kpfi = — kz, and (r, ip) form a system of cylindri- 
cal coordinates attached to the rod axis. 

The field components along z can be expanded into 



cylindrical harmonics. For example, 

oo 

E anJin\ {kp,mr) eJ•"^e-^■'=-^ r<R 
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E KH\^{ {kp^or) eJ"'^e-^ '=^^r > R 

n= — oo 

(2) 

where a„ and 6„ are the unknown coefficients of the ex- 



pansion. H; 



(2) 



Jn —jYn is the Hankel function of second 



kind and order n, and kp^m = —j^/^z — /3^£. The field 
Hz has a similar expansion. 

The transverse fields E|[ and H|| (projections into the 
xoy plane) can be written in terms of Ez and Hz'. 
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{-jkz'^llEz+jPuz X ViWoHz) (3) 



£ and kp,i = kp^m inside the rod, and 



?7oH|| = -2- {-jPsiUz X V\\Ez-jkz ^wVoHz) (4) 

In the above, % is the free-space impedance, V|| = 
^ Q\ i 

£i = 1 and kp^i = kp^ outside. 

The unknown coefficients can be calculated by impos- 
ing the continuity of the tangential electromagnetic fields 
(i.e. Ez, E^, Hz, H^) at the interface r = R. 

The {z component of the) electric dipole moment (per 
unit of length) is given by. 
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and so the polarizability per unit of length is. 



= (e - 1) 27ri? 



x=y=0 

ao Ji (kp.mR) 



En 



(6) 



This result shows that ao is the unique coefficient re- 
quired to calculate the polarizability. Since Maxwell 
equations are separable in cylindrical coordinates, to cal- 
culate ao it is sufficient to impose the boundary condi- 
tions to the terms associated with the cylindrical har- 
monic n = 0. It can be easily verified that the n = 
term of the magnetic field Hz vanishes (this happens be- 
cause iJ™"^ = and ^ = for the n = harmonic). 
Thus, only the tangential components Ez and H^ have 
non-trivial n = coefficients. Using (2)-(2) in (3)-(4), 
and imposing the continuity of the tangential fields, we 
find that, 

ao Jo {kp,mR) = hoH^o^ [kp^R) + EqJo {kp,oR) (7) 
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where the prime " "' denotes the derivative of a function. 
Solving for ag wc readily obtain: 
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Substituting this result in (6), we obtain the desired 
electric polarizability: 
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(10) 



In this work, we consider that the radius of the rods 
is always much smaller than the wavelength, or equiva- 
lently, Rkp^ « 1. In these circumstances, (10) simpli- 
fies to, 
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where C is the Euler constant. It is important to note 
that because the rods are infinitely long, the polarizabil- 
ity depends not only on the frequency of the incoming 
wave, but also on the wave vector component kz- If e 
approaches —oo the above result reduces to the case of 
perfectly conducting rods studied in [14]. 

Proceeding similarly, we can obtain the well-known for- 
mula for the electric polarizability (per unit of length) in 
the transverse plane {xoy plane): 
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(12) 



Provided the permittivity of the rods is not too close to 
the resonance e = — 1 and the rods are very thin, the 
polarizability can be neglected in the transverse plane. 



III. EFFECTIVE PERMITTIVITY MODEL 

In the following, we derive an effective permittivity 
model for the medium formed by a periodic array of ENG 
rods. As depicted in Fig. 1, the rods are arranged in a 
square lattice and the spacing between the rods (lattice 
constant) is a. As is well-known, each electromagnetic 
(Floquet) mode in a periodic medium can be associated 
with a wave vector k = {k^, ky,kz). For convenience, we 
define k|| = {kx,ky,Qi). 

To compute the effective permittivity we use the mix- 
ing formula : 



= = 1 = 
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(13) 



where Aceii = a?, I is the identity dyadic, the super- 
script "-1" represents the inverse dyadic, and Cint is the 
interaction dyadic calculated in Appendix B. The mix- 
ing formula is derived in Appendix A, and is valid under 
the condition that the dimensions of the cross-section 
are much smaller than the lattice constant, and that 
|k|||a << 2n. As discussed in Appendix A, even though 
(13) reminds Clausius-Mossotti formula, things arc not 
so plain because the lattice has some intrinsic dispersion. 
To keep the readability of the paper the details have been 
moved to Appendix A. 

Next we substitute (11) and (12) in (13). Using (B4), 
we find that the effective permittivity in the transverse 
plane is 
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where fv = nR^/Acdi is the volume fraction of the rods. 
Provided the permittivity of the rods does not satisfy e « 
— 1 and the rods are very thin, we can assume that e^x = 
Cyy ~ 1. For simplicity, we shall assume this situation in 
the rest of the paper. On the other hand, from (13) it is 
clear that, 



= 1 + 
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where Cint,zz is given by (B7). So, after further simplifi- 
cations, we obtain that: 



£zz = 1 



1 _ 

(e-l)/v 



(16) 



where (3p is the plasma wave number defined consistently 
with the results of [14] for perfectly conducting wires: 
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Note that in general e is a function of frequency. Formula 
(16) gives the effective permittivity of an array of diluted 
ENG rods. The first important observation is that the 
medium is spatially dispersive. Indeed the permittivity 
depends not only on the frequency fi = oj / c, but also on 
the component of the wave vector parallel to the rods. 
This means that the medium is nonlocal, i.e. in the spa- 
tial domain the electric displacement vector and the elec- 
tric field are related through a spatial convolution rather 
than by a multiplication. Secondly, we note that if the 
permittivity of the rods approaches —00, (16) reduces to 
the formula derived in [10] for the wire medium, consis- 
tently with the observation made in the introduction of 
this paper. Finally, if we put /? = and assume on-plane 
propagation, i.e. fcz = 0, the effective permittivity sim- 
plifies to Ezz = 1 + (e — 1) /y, which is the exact formula 
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in the static the Hmit for non-dispersive dielectrics [16]. 
Thus, very interestingly, formula (16) is in a certain sense 
the average of these two limit situations. Note that even 
though in this work our main interest is the analysis of 
ENG rods, the proposed model is also valid for dielectrics 
with positive real part of the permittivity. In the next 
sections, we characterize the electromagnetic modes sup- 
ported by the rod medium and validate the model with 
numerical simulations. 





IV. CHARACTERIZATION OF THE 
ELECTROMAGNETIC MODES 



FIG. 2: Reflection of a plane wave by a semi-infinite rod 
medium. a)The axes of the rods are parallel to the interface, 
b) The axes of the rods are normal to the interface. 



It is evident that the waves in the homogenized 
medium can be decomposed into transverse electric (TE) 
modes (electric field is normal to the axes of the ENG 
rods), and transverse magnetic (TM) modes (magnetic 
field is normal to the axes of the ENG rods). As ex- 
plained in the previous section, we shall assume in this 
paper that exx = £yy ~ 1 (i-e. that the rods are very thin, 
and that the permittivity of rods is not close to e « — 1). 
Within this approximation, the TE-modes do not inter- 
act with the rods, and thus their dispersion characteristic 
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and the associated average electric field is (the propaga- 
tion factor e^^^ '' is implicit) 



SI 



(19) 



On the other hand, the TM-modes satisfy the character- 
istic equation: 



iP' - ki) 



(20) 



The above equation cannot be solved explicitly as a func- 
tion of (3, because the permittivity of the rods is itself a 
function of [3. The corresponding average electric field is 
(for fc, ^ 0): 
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(21) 



The associated magnetic field can be calculated using 
(A2). 

To understand better the nature of the TM-modes, 
next we study a reflection problem. Let us consider a 
semi-infinite rod medium illuminated from the air side 
with a plane wave. We analyze two different geometries, 
as depicted in Fig. 2. 

Firstly, let us suppose that the axes of the rods are par- 
allel to the interface x = (Fig. 2a). The incident wave 

vector is k™'^ = {—jlo, ky,kz) with 70 = y/fc^ + — (3'^. 

It is well-known that the component of the incident wave 
vector parallel to the interface, (0, fcjy,/;^), is preserved. 



In this case only one TM-mode is excited in the artificial 
medium (besides the TE-mode). Indeed, from (20) the 
component kx of the wave vector inside the rod medium 
is given by: 



rod 



(22) 



Notice that the right-hand side of the above equation 
only depends on the geometry and parameters of the 
medium, on the wave number P of the incident wave, 
and on the components of the incident wave vector that 
are preserved (0, fcy, fc^). 

Next suppose that the rods arc normal to the interface 
z — (Fig. 2b). The incident wave vector is now k™"^ = 

{kx,ky, ~j^o) with 7o = ^Jk^~+k^~~P^ . The interesting 

thing is that for this configuration two TM-modes can 
be excited inside the rod medium. Indeed, since k|| = 
{kx,ky,0) to find the excited electromagnetic modes one 



needs to solve (20) for kz 
using (16), show that: 



kl^ (3^ 



Straightforward calculations. 
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where we defined the parameter (3c = Pd^) as. 
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(24) 



Note that provided the permittivity of the rods is less 
than the permittivity of the host medium, /3c is a positive 
real number (with the same unities as /3; for simplicity 
the rods are assumed lossless, otherwise /?c becomes a 
complex number). Also, (3c is in general a function of 
frequency since e also is. From (23) it is seen that there 
are two different solutions for fcz, and hence two TM- 
modes, besides the TE-mode, can propagate inside the 
artificial medium. This phenomenon is a manifestation 
of spatial dispersion, and is also characteristic of the wire 
medium [10]. We also note that the average electric field 
for both TM-modes is calculated using (21). 
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FIG. 3: Plot of kz as a function of normalized frequency /3a, 
for R = 0.05a and fe|| = 0. The permittivity e follows a 
Drude type model (see the text). The solid line represents 
the "exact" result, whereas the dashed line represents the 
results calculated with the proposed model. 



To illustrate the discussion, we plot in Fig. 3, kz as 
a function of normalized frequency /3a for the parame- 
ters R — 0.05a, =0, assuming that the permittivity 
follows the Drude model e = 1 — /3^//3^ with (normal- 
ized) plasma wave number /3ma = 12.0. The dashed line 
curve represents the results calculated using (23). The 
solid hue curve corresponds to the data calculated by 
substituting (10) in (15) (with no approximations) and 
calculating Cn-it,zz using (B5) with r = r' = 0, but with- 
out making any assumptions with respect to k||a being 
negligible. Hence, the permittivity becomes a function 
of not only (3 and kz, but also of the other components 
of the wave vector. The permittivity function obtained 
in this way is substituted in (20), and the corresponding 
equation is solved numerically. A similar procedure was 
used in [9, 14], and so further details are omitted here. 
These results can be regarded as "exact" within the thin 
rod approximation. As seen in Fig. 3, the agreement 
for the first TM-mode is always excellent. On the other 
hand, the second TM-mode is not so accurately predicted 
for relatively small wavelengths. The reason is twofold. 
The first reason is that for larger frequencies the long 
wavelength limit approximation is not so accurate. The 
second reason will be discussed later. Fig. 3 shows that 
for = one of the electromagnetic modes propagates 
with the speed of light. It is also important to refer that 
for very low frequencies one of the TM-modes is cut- 
off (complex imaginary propagation constant). This is 
because the ENG rods behave as perfect conductors in 
the static limit (when the permittivity follows the Drude 
model) . 

To give further insight about the TM-modes, let us 
study different limit situations. First, suppose that at 
some frequency (3c{ijj) « /3p, i.e. the permittivity of 
the rods is very large in absolute value. In this case, 
(23) reveals that one of the modes has the dispersion 
characteristic k1 = 0^ , and that the other mode has the 



dispersion k1 = 0^ — 0^ — k'^y The former mode can be 
readily identified with the well-known transverse electro- 
magnetic (TEM) dispersionless mode of the wire medium 
(perfectly conducting wires), while the latter is the TM- 
mode of the wire medium. 

Consider now the case k|| « 0, i.e. paraxial incidence. 
Using a Taylor expansion we obtain: 
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(25) 



The above formula shows that if either pc{(^) << Pp or 
/3c{w) >> /3p, one of the modes becomes dispersionless 
with respect to k|| (i.e. the coefhcient associated with 
k^ vanishes). The former case was already discussed. 
As to the latter case, the pertinent mode has dispersion 
k^ « P'^ + 01. But this implies that kz > P and thus 
this mode is a surface wave guided along the rods. In 
[17] it was proved that a ENG rod is able to support 
tightly bounded surface modes that propagate electro- 
magnetic energy with subwavelength beam radius. For 
metallic materials the surface modes are surface plasmon 
polaritons. It was shown that the energy becomes more 
confined to the vicinity of the dielectric waveguide when 
the effective index of refraction n^ff — kz/P increases. 
This important result justifies why in the rod medium 
one of the TM-modes becomes independent of k||. In 
fact, when (3c{(^) >> Pp each guided mode is confined 
to a small vicinity of the respective ENG rod, there is 
no interaction or coupling between the rods, and conse- 
quently one of the TM-modes becomes dispersionless. As 
referred before, when Pd^^) << Pp there is also a quasi- 
TEM dispersionless mode. However this mode is qual- 
itatively very different from the mode that arises when 
/?c(w) >> Pp. Indeed, while in the latter case the en- 
ergy is propagated tightly bounded to the ENG rods, in 
the TEM mode case the field energy is distributed more 
or less uniformly by the whole volume. As discussed in 
the introduction, the dispersionless modes may be used 
to canalize the electromagnetic radiation through the rod 
medium and achieve sub-wavelength imaging. A detailed 
analysis of this topic is out of the scope of the present 
paper, and will be reported elsewhere. We also refer that 
the other TM-mode, still assuming that Pc{oj) » Pp, 
has to a first approximation the dispersion characteristic 
kz~P'^ — k'^y i.e. approximately the same dispersion as 
the TE-mode. 

In Figs. 4 and 5, kz is plotted as a function of k^ for 
R = 0.05a and several different values of the normalized 
frequency Pa (the permittivity of the rods follows the 
same Drude model as before). For convenience we show 
the results in two different figures. For very long wave- 
lengths /3c (^) << Pp, and consequently one of the modes 
is cut-off. The dispersion of the other mode (shown in 
Fig. 4) becomes increasingly flat as the frequency (and 
consequently /3c) decreases. Note that from (24) and as- 
suming a Drude type model, /3c increases monotonically 
with the frequency. 
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FIG. 4: Plot of fcz as a function of k^a for R — 0.05a and ky — 
0. Contours /3a = 0.1, 0.3, 0.5, 0.7, 0.9 for the TM-modes. The 
permittivity e follows a Drude type model (see the text). 




FIG. 5: Plot of fcz as a function of kxCi for R = 0.05a and 
ky = 0. Contours /3a = 1.0, 1.2, 1.4, 1.6, 1.8, 2.0 for the TM- 
modes. The permittivity e follows a Drude type model (see 
the text). 



At some point, as the frequency increases, the TM- 
mode that was cut-off starts to propagate. In this case, 
(Fig. 5), for each fixed frequency there are two differ- 
ent contours, i.e. two propagating TM-modes. As the 
frequency increases even more, /3c becomes comparable 
or larger than Pp, and the band structure of one of the 
TM-modes becomes practically flat, consistently with the 
fact that the energy propagated by this mode is tightly 
confined to the vicinity of the ENG rods. 

In Figs. 4 and 5 it can also be seen that around fc|| « 
the wave normal contours of one of the modes are to 
some extent hyperbolic. In fact, using (25) it can be 
easily checked that one of the modes has dispersion char- 
acteristic such that the signs of the coefficients associated 
with and fc^ are symmetric. This property is more im- 
portant for frequencies such that /3c(w) ~ Pp. Note that 
if the medium was anisotropic, with no spatial disper- 
sion, and negative permittivity along the z-direction the 
contours would also be hyperbolic. It is well-known that 



hyperbolic contours may originate negative refraction at 
an interface. 

It is also important to refer that in the limit e — > 1, 
the ratio Pd^)/ Pp becomes infinitely large (see (24)) and 
consequently kz/ P is also infinitely large. This means 
that in these circumstances the localized TM-mode either 
cannot be excited by an incoming wave, or if it is excited 
it is killed by losses. Thus, only the other TM-mode, 
with dispersion k1 k, /3^ — fc^, will propagate. This is 
consistent with the fact that as e — > 1 the medium shall 
have the same properties as free-space. 

To conclude this section, we will discuss the scope of 
application of the permittivity model (16). We remember 
that the results were derived for thin rods, R << a, and 
under the assumption that |A:p.o|a << tt and Pa << n. In 
general the modes that propagate in the long wavelength 
limit satisfy the previous conditions without problems. 
However, there is one exception at the problem at hand. 
In fact, when Pci^jj) » Pp the radial constant kp^ be- 
comes complex imaginary for the TM-mode associated 
with the surface mode (surface plasmon) . As discussed in 
[17], when the beam radius is subwavelength the effective 
index of refraction rieff = k^/P becomes very large, and 
in that case the condition |A:p.o|a << tt may not be ob- 
served. This situation affects the accuracy of our model 
when Pc{(jj) » Pp- Indeed, the error in the TM-mode 
associated with the surface mode becomes non-negligible 
in this situation. The other TM-mode is still accurately 
predicted. This result justifies deterioration of the agree- 
ment in Fig. 3, as the frequency (and consequently, for 
the Drude Model, also Pc) increases. 

Fortunately, it is easy to solve this problem. In fact, 
when Pc{i^) » Pp the dispersion characteristic of the 
pertinent TM-mode is essentially the same as the disper- 
sion characteristic of the guided mode supported by a 
single ENG rod. This dispersion characteristic is deter- 
mined in [17], and is equivalent to the condition — 0, 
where a~} is given by (10). Thus, to summarize our find- 
ings, the TM-modes can be accurately calculated using 
(23), except when Pc{(^) » Pp which yields less accurate 
results for the mode with higher kz ■ In this case the cor- 
responding TM-mode is dispersionless, and follows the 
same characteristic as the guided mode supported by a 
single dielectric rod. 



V. THE REFLECTION PROBLEM WITH RODS 
PARALLEL TO THE INTERFACES 

To further validate the proposed permittivity model, 
we will study the reflection of electromagnetic waves by 
a rod medium slab with finite thickness. We will suppose 
that the rods are parallel to the interface (see Fig. 6). 
The slab consists of Nl layers of rods. The structure is 
periodic in y and z, and the dielectric rods stand in free- 
space. The interfaces x = x'j^ and x = x'p^ = x'j^ + N^a 
are represented by the dashed lines (since the rods stand 
in air the definition of the interfaces is a bit ambiguous; 
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this will be discussed below with more detail). The rods 
in the leftmost layer are in the plane x = Xq. 

As discussed in the previous section, when the rods are 
parallel to the interface, equation (22) has only one solu- 
tion for k^. For simplicity, we will restrict our attention 
to the case in which either — (Fig. 6a) or ky — 
(Fig. 6b). For these particular geometries, an incident 
plane wave polarized as depicted in Fig. 6 can only excite 
the TM-mode inside the rod medium. As is well-known 
[18], in the general case where ky and k^ are simultane- 
ously different from zero, both the TM- and TE-modes 
are excited (the medium is birefringent). Apart from the 
more heavy notation, the general case poses no additional 
difficulties. 
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FIG. 6: A plane wave illuminates a slab of the rod medium 
{Nl = 2). The rods are parallel to the plane x = 0. a) fcz = 0. 
b) ky = 0. 



of the reflection coefficient and for the transmission co- 
efficient. The previous results also demonstrate that the 
homogenization model is useful to study not only incident 
propagating plane waves, but also part of the evanescent 
spectrum. 
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FIG. 7: Reflection coefficient as a function of the wave vector 
component parallel to the interface. The slab consists of Nl = 
5 layers, and the rods have R = 0.05a, and e = — 30 at /3a = 
1.0 The solid line represents the full wave MoM data, and the 
dashed line represents the data obtained using the analytical 
model. 



Using (A2), (21), and (22) and matching the tangen- 
tial components of the electric and magnetic fields at the 
interfaces, we find that the reflection coefficient referred 
to the plane x = x'^ is given by: 



tanh(7™rf) (to -7,^ 



P = 



27o7m + tanh (7„d) (7^ + 7,^) 
tanh (jrnd) 



(e'7o 



2£7o7m + tanh (7„d) (£2-^2 _(_ 



(26) 
(27) 



where d = Ni^a is the thickness of the slab, 7™ = j k^ rodi 
k1 — /?2, and (26) corresponds to Fig. 6a) 



70 



1,2 

Ky 



with kz — 0, and (27) corresponds to Fig. 6b) with 
ky = 0. 

Next, in order to demonstrate the accuracy of the the- 
oretical results, the analytical model is tested against full 
wave data computed with the periodic moment method 
(MoM) [19]. In the first example we consider that 
R = 0.05a, and e = —30.0 at (3a — 1.0 (for simplicity, 
losses are neglected) . A plane wave polarized as depicted 
in Fig. 6a illuminates 5 layers of rods {Nl = 5). The am- 
plitude of the reflection coefficient is depicted in Fig. 7 
as a function of kya. Note that the angle of incidence is 
such that sin(^ = kyjfj. For ky > (3 the incident wave is 
evanescent. The solid line represents the MoM full wave 
data. The dashed line represents the results computed 
using the proposed permittivity model (formula (26)). It 
is seen that that the agreement between the two sets of 
data is good. Similar agreement is obtained for the phase 
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FIG. 8: Reflection coefficient as a function of kya for Nl ~ 1 
layer, R = 0.05a, and e = —30 at /?a = 1.0. The solid and 
dashed lines are defined as in Fig. 7. 

As noted before, since the rods stand in free-space the 
position of the interfaces and thickness of the slab are a 
bit ambiguous. Notice that the thickness of the homog- 
enized slab was taken equal d = N^a, apparently with 
good results. Next, to test if this choice still yields accu- 
rate results for very thin slabs, the reflection coefficient 
is computed for the same structure, except that now the 
slab has only one layer of rods {N^ — 1). The refiec- 
tion coefficient is depicted in Fig. 8 and has a peak at 
kya — 1.0, which corresponds to the transition between 
propagating waves and evanescent waves. As seen, even 
though the slab is so thin the agreement is still remark- 
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ably good. This is a bit a surprising, because for such a 
thin slabs one would expect that the interface effects and 
granularity of the artificial medium would prohibit the 
homogenization of the structure using the bulk medium 
average fields. 

In the next example, we study what happens if the an- 
gle of incidence Lp is kept constant {if = 45[deg]), and the 
frequency is varied. Now R = O.Ola, Nl = 5, and the 
permittivity follows the Drude model e ~ 1 — /3^//3^ with 
plasma wave number PmO, = 12.0 or /^mO = 80.0. The 
calculated results are shown in Fig. 9. For relatively low 
frequencies the two sets of data agree very well, but as 
the frequency increases the agreement progressively de- 
teriorates, since the long wavelength limit approximation 
is no longer valid. Notice that for relatively low frequen- 
cies the medium blocks the incident radiation because 
the rods effectively behave as perfectly conducting wires. 




12 3 4 

Normalized Wave Number, pa 

FIG. 9: Reflection coefficient as a function of the free-space 
wave number (3a for R = O.Ola, and Nl ~ 5 layers. The 
permittivity of the rods follows a Drude-type model (see the 
text). The solid and dashed lines are deflned as in Fig. 7. 
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FIG. 10: Reflection coefflcient as a function of kza for Nl ~ 5 
layer, R = 0.05a, and e = —30 at /3a = 1.0. The solid and 
dashed lines are defined as in Fig. 7. 

So far the wave vector of the incident wave was al- 
ways perpendicular to the axes of the ENG rods, and so 



the effects of spatial dispersion were hidden. In the last 
example, we consider the very different propagation sce- 
nario depicted in Fig. 6b. The parameters of the rods 
are R — 0.05a, e — —30, and Nl — 5. The reflection 
coefficient for f3a — 1.0 is shown in Fig. 10 as a function 
of k^a. Note that the angle of incidence 6 for propagat- 
ing waves satisfies sin 9 = kz/(3. Fig. 10 shows that the 
agreement between the numerical and analytical results 
is still quite satisfactory, except near the transition be- 
tween the propagating waves and the evanescent waves 
(k.a^l.O). 

It is important to refer that the configurations stud- 
ied in this section (rods parallel to the interface) are not 
appropriate to achieve sub- wavelength imaging. For that 
application the rods must be perpendicular to the inter- 
face, as shown in Fig. 2b. As discussed in section IV, 
for this geometry three modes are excited in the (homog- 
enized) rod medium, and thus an additional boundary 
condition is necessary to solve the scattering problem. 
The same situation occurs in the wire medium [15]. The 
analysis of this problem is out of the scope of this work 
and will be presented elsewhere. 



VI. CONCLUSION 

We studied the electrodynamics of a periodic array of 
thin ENG rods. Using the polarizability of a single rod 
and integral equation methods, we derived new nonlo- 
cal permittivity model for the artificial medium that ac- 
curately describes the propagation of waves in the long 
wavelength limit. We discussed the effects of spatial dis- 
persion in the context of the reflection problem. It was 
proved that when the rods are parallel to interface only 
two modes (TE and TM) can be excited in the artifi- 
cial medium. However, when the axes of the rods are 
normal to the interface, two TM-modes, besides the TE- 
mode, can be excited in the medium, as a manifestation 
of spatial dispersion. It was demonstrated that the wave 
normal contours of one of the TM-modes are intrinsically 
hyperbolic. It was proved that the rod medium supports 
dispersionless modes that propagate along the axes of the 
rods, and it was speculated that this important prop- 
erty may allow sub-wavelength imaging of electromag- 
netic waves at the infrared and optical domains, using 
the idea proposed in [4]. It was shown that the energy 
of the dispersionless modes can be loosely bounded to 
the ENG rods (in this case the wave is essentially trans- 
verse electromagnetic), or alternatively tightly confined 
to a vicinity of the rods. In the latter case, the modes 
have approximately the same dispersion as the guided 
surface mode supported by an individual rod. The re- 
flection problem was investigated in detail for the case 
where rods are parallel to the interfaces. The developed 
theory was successively tested against full wave data cal- 
culated with the MoM. 
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In the above, $p = $p(r|r') is the dynamic potential 
created by a phase-shifted array of line sources, 

VHp + pHp = - ^ ^ (r|| - rj| - ri)e-^''- ('-•') (A8) 



APPENDIX A: DERIVATION OF THE MIXING 
FORMULA 

In this Appendix we derive the mixing formula (13) 
used to homogenize the rod medium. 

Let us consider a generic electromagnetic Floquet 
mode (E,H), associated with the wave vector k = 
{kx, ky,kz) and the wave number = ui/c, i.e. the fields 
satisfy Maxwell-Equations and (E, H) exp (jk.r) is a pe- 
riodic function in the lattice. The average electric field 
is defined as, 

Eav = y E (r) e+^" '^ '^d^'r (Al) 



and Hav is defined similarly. In the above, Q represents 
the unit cell of the periodic medium. Prom Maxwell- 
Equations it can be proved that: 



-k X Eav + P %Hav = 



/? ( Eav + — ) + k X %Hav = 



(A2) 
(A3) 



where the generalized polarization vector is given by, 

P 

£0 



= / (£ - 1) E ^ "d^r (A4) 

Vcell J 



Straightforward calculations show that these relations 
imply that: 



_ - kk P 



(A5) 



The above equations are exact and completely general. 

Next we will apply these results to the rod medium under 
study. To begin with, we note that since the crystal is 
invariant to translations along the z-direction the fields 
depend on z as exp{—jkzz). Using standard Green func- 
tion methods [18], it can be proved that the electric field 
has the following integral representation, 



E(r) = (£-!) J /32Gp(r|r').E(r') dV 
s 



(A6) 



where the primed and unprimed coordinates represent 

the source and observation points, respectively, S = 
{ (x', y',0) : x'^ -I- y'^ < R^} is the cross-section of the di- 
electric rod in the unit cell, and the Green function dyadic 
is defined by. 



Gp=\^I+ — (V|| - jkzUz) (V|| - jkzUz) ] % (A7) 



where S is the Dirac function, I = (ii,i2) is a double 
index of integers, ri = a(ii,i2,0) is a lattice point, and 
rj| = {x,y,0). Thus $p is intrinsically two-dimensional, 

depending on z and z' as e Furthermore, it is 

obvious that the Green function only depends on the rel- 
ative position r|| — r||. Note that the Green function can 
be written as a superimposition of the potentials created 
by the line sources: 



(A9) 



where $0 



H^'> (/cp,o|r|| -r'|||)/4j is the 



potential created by a line source placed at r'l, i.e. the 
solution of (A8) when the summation in the right-hand 
side is restricted to the index 1 = 0. Physically, (A6) 
establishes that the electric field at some point of space 
is the superimposition of the fields radiated by all the 
dielectric rods of the lattice. 

The Green potential can also be written as a Fourier 
series since it is a (pseudo) periodic function of the wave 
vector: 



g-jkj.(r-r') 



1 g— J'^J•^,'— ' J 



(AlO) 



where Aceii = a^, J = (ji, j2) is a double index of inte- 
gers, kj = k-l- k5, and kj = 27r/a(ji, j2,0). 

Now that the necessary theoretical formalism was in- 
troduced, we are ready to study the homogenization 
problem in the rod medium. To begin with, we note that 
from (A9) and (AlO) the Green potential is singular in 
the spatial domain, i.e. when r|| — rj| « (source region), 

as well as in the spectral domain, i.e. when k.k « 
(long wavelength limit). Since the integral (A6) is de- 
fined over the source region and we want to study the 
electromagnetic modes that propagate in the long wave- 
length limit, it is convenient to single out the terms that 
make the Green function singular and decompose it as 
follows: 



$p = $0 + 



-jk.(r-r') 



Acell k.k-/?2 



reg 



(All) 



where ^reg, which is defined implicitly by the above for- 
mula, is a regular function both in the spatial domain 



(source region r 



Ml 



0) and in the spectral do- 



main (long wavelength limit). Using this decomposition 
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in (A6) we find that: 

E (r) = E.ve-^''-'- + " 1) / Z^'^o ( r| r') .E d^r' 

s 

+ i^-'^)J P^^reg {r\r').E(fr' (A12) 

s 

where Eav is the average field in the crystal, and Gq and 
Greg are defined as Gp, except that the Green potential 
is replaced by <E>o and ^reg, respectively. To obtain the 
above result we used (A5). Next, we use the fact that the 
ENG rods are assumed to be very thin (R/a << 1), and 
that we want to investigate the electrodynamics of modes 
that propagate in the long wavelength limit (|k|| |a << 27r 

and 13a « 2tt). Since the dyadic Greg is regular in both 
the spatial and spectral domains, it is legitim to write 
(putting r|| = rJi = and k|| = {kx,ky,0) = 0), 

^reg (r| r'; k, 13) « ^reg {z\z' ; k„ (3) (A13) 

being the formula valid in the source region. For conve- 
nience, we introduce the following interaction dyadic: 

^nt=/3'Sreg(0;fc„/?) (A14) 

Then, it is clear from (A12) that, 

E(r)« (^Eav+ant.|^)e-^'=^^ 

+ {e-l) J l3^%{r\r').Ed^r' (A15) 
s 

provided the observation point r is near the source region 

and |k|||a << 27r. In above, we introduced the electric 
dipole moment (per unit of length), p, of the dielectric 
rod in the unit cell. Now, the key result is that (A15) 
is formally equivalent to the (integral) equation obtained 
when a single rod is illuminated by a plane wave with 

electric field with amplitude = ^Eav + Cint-^^ and 

the same wave vector component along the z-direction. 
In other words, when the rod in the unit cell stands alone 
in free-space and is illuminated with the above defined 
plane wave, the total electric field also satisfies (to a first 
approximation) (A15) in the source region. But this re- 
markable result implies that: 

^ =Se. fEav+Oinf-) (A16) 

where Txe is the electric polarizability tensor for a single 
rod. The term inside brackets in the right-hand side can 
be regarded as the local field that polarizes a single rod 
embedded in the dielectric crystal. Within the thin rod 
condition and for long wavelengths, the above solution is 
exact. 



We are now ready to calculate the effective permit- 
tivity dyadic. Since the (macroscopic) electric displace- 
ment vector D is given by D = £oEav + p/Aceii, and 
the effective permittivity tensor must guarantee that 
D = eoS-Eav, we conclude that the effective permittivity 
of the rod medium is given by the mixing formula (13). 
Note that (13) is completely general and is valid indepen- 
dently of the specific geometry of the transverse section 
of the rod. 

At this point it is appropriate to compare (13) with 
the classic homogenization approach. It is striking that 
(13) reminds Clausius-Mossotti formula [16, 18]. In- 
deed, if we could identify the interaction dyadic Cint with 
l/2Aceii the formulas would be the same (note also that 
the rods are arranged in a square lattice). In Appendix 
B we calculate the interaction dyadic in closed-form us- 
ing the static limit approximation. Equation (B4) shows 
that the interaction dyadic is different from 1 /2Accii only 
along the z-direction. This important result shows that 
Clausius-Mossotti formula is not valid for media with 
cylindrical inclusions. More specifically it fails to predict 
the effective permittivity along the direction in which the 
crystal is uniform. This is an indirect manifestation of 
spatial dispersion. 

We also mention that the interaction dyadic defined 
by (A14) is not equivalent to the dynamic interaction 
constant defined in other works (see for example [20]). 
Indeed, in our definition we extracted the singularities 
in both the spatial and spectral domains, while other 
works usually only extract the singularity in the spatial 
domain. It is clear from our previous analysis that it 
is the singularity in the spectral domain that indirectly 
defines the relation between the local field that polarizes 
the rod and the average field. 

APPENDIX B: CALCULATION OF THE 
INTERACTION DYADIC 

Here we calculate the interaction dyadic defined by 
(A14). It can be written as: 

Cint = (0+ (V|| -jfc^U,) (V|| -iA;,U,)) $reg (Bl) 

where the right-hand side of the expression is evaluated 
at r = r' = and k|| = 0. From (A8) and (All) it is 

clear that: 



V^^reg + P'^^reg = 




Putting r = r' = and k = in the above equation, and 
letting (3 approach zero (static limit), we find that: 

V^$reg(r = r' = 0;k = 0)| =J- (B3) 
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Moreover, because of the symmetry of the square lattice 

and k = the following 

= 0, and 



it is evident that if r 

dxidx 



Using (A9) and (All), and putting r = r' = and k|| = 
0, we obtain that: 



relations hold, 



if « J, -g^ 



So using (B3), we conclude that in the Cint,zz = {P^ - k^) ^ —H^^^ {kp^ |ri|) + j— (B6) 
lit fk = and 8 = 0) the interaction dvadic is 



static limit (k = and /? = 0) the interaction dyadic is 
given by: 



(B4) 



The above result is consistent in the xoy plane with the 
(two dimensional version of the) Clausius-Mossotti for- 
mula. However, along the z-dircction, maybe a bit sur- 
prisingly, the interaction constant vanishes in the static 
limit. Next, we will estimate Cint,zz in the dynamic case. 
From (Bl), we have that: 



The series in the right-hand side was evaluated in [14]. 
Using the results of [14], and assuming that kpfiU « it, 
we obtain: 



int.- ~ '^P.O U ■ 277 - V 4^ y ' 2tT 
1 1 



12 ^ 7r|n| e2'^l"l - 1 

n—l ' ' , 



(B7) 



a. 



(B5) 



where C is the Euler constant. 
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We propose an effective permittivity modei to homogenize an array of long tliin e-negative rods 
arranged in a periodic lattice, ft is proved that the effect of spatial dispersion in this electromagnetic 
crystal cannot be neglected, and that the medium supports dispersionless modes that guide the 
energy along the rod axes. It is suggested that this effect may be used to achieve sub-wavelength 
imaging at the infrared and optical domains. The reflection problem is studied in detail for the case 
in which the rods are parallel to the interfaces. Full wave numerical simulations demonstrate the 
validity and accuracy of the new model. 
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I. INTRODUCTION 

In recent years there has been a lot of interest in the 
propagation of electromagnetic waves in artificial mate- 
rials, and particularly in materials with negative index of 
refraction [? ? ]. The research has been mainly driven 
by the possibility of using these materials to miniatur- 
ize several devices and waveguides, and develop "perfect 
lenses" able of focusing electromagnetic radiation with 
sub- wavelength resolution [? ]. Recently, a different 
approach to achieve sub-wavelength resolution was pro- 
posed in [? ] and demonstrated experimentally in [? ]. 
In [? ] the idea is to use an artificial material formed by 
an array of perfectly conducting wires (wire medium) to 
guide the waves from the input plane to the output plane, 
and then reconstruct the image "pixel by pixel", explor- 
ing a Fabry-Perot resonance that in the wire medium 
occurs simultaneously for all the spatial harmonics. The 
resolution of this transmission device is only limited by 
the lattice constant, i.e. by the spacing between the 
wires. The problem with the configuration studied in [? 
] is that it is limited to the microwave regime because in 
the optical domain perfect electric conducting materials 
are not available. Nevertheless, we will suggest in this pa- 
per that it may still be possible to use the same concept to 
achieve sub- wavelength resolution at the infrared and op- 
tical domains, provided e-negative (ENG) rods are used 
to guide the electromagnetic radiation instead of metallic 
wires. At the infrared and optical frequencies all metallic 
materials have permittivity with negative real part. It is 
known that their dielectric constant is well represented 
by the Drudc model e = 1 — iJ^juj^ , where lo^ is the 
plasma frequency (for simplicity the lossless model was 
considered). Thus, we envision that either silver or gold 
rods may used to fabricate transmission devices that are 
able to propagate the subwavelength information of an 
image. 

With this motivation, we will investigate in this paper 
the propagation of electromagnetic waves in an artificial 
medium formed by thin ENG rods (see Fig. 1), and we 
will show that the structure can be homogenized and de- 
scribed by an effective permittivity tensor provided spa- 



tial dispersion is taken into account. We will show that 
our homogenization model predicts that the rod medium 
may support a mode that propagates along the axes of 
the rods with the same phase velocity, independently of 
the transverse wave vector. As proved in [? ], this is a 
key property to operate the material in the canalization 
regime and achieve sub-wavelength resolution. 




FIG. 1: Medium formed by long thin ENG rods arranged in 
a square lattice. 

Previous works related with the homogenization of 
ENG rods are discussed next. In [? ] the oblique 
propagation of electromagnetic waves through an array 
of aligned fibres was examined, and the associated prob- 
lem of numerical instability discussed. In [? ] a method 
was proposed to compute analytically the band structure 
of wire mesh crystals. In [? ] a model was proposed to 
homogenize a medium with metallic rods, but only the 
on-plane case was considered. In [? ] an homogenization 
model was proposed to characterize a structure similar 
to the one studied in this paper. It was demonstrated 
that the artificial material was characterized by spatial 
dispersion, and that the band structure of the photonic 
crystal had several branches. However, the results of [? 
] are restricted to the case in which the permittivity of 
the rods follows a Drude-type plasma model, and besides 
that the derived formulas for the effective permittivity 
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arc very cumbersome, lead to non-analytical dispersion 
characteristics, and more importantly hide the physics 
of the problem. In this paper, we will derive a simpler 
and more intuitive model, describe new phenomena, and 
prove that the new approach can characterize accurately 
the electrodynamics of the artificial medium. 

The homogenization of the wire medium is also closely 
related with the subject under study. In fact, the case of 
perfectly conducting wires can be regarded as the limit 
situation in which the permittivity of the ENG rods is 
—00. In [? ] an homogenization model was derived for 
the wire medium, and it was proved that this artificial 
material suffers from strong spatial dispersion even for 
very long wavelengths. Later, in [? ? ], these results 
were generalized for 2D- and 3D lattices of connected 
and unconnected wires. In [? ? ] the reflection problem 
in a wire medium slab was investigated, and in [? ] it was 
proved that in general an additional boundary condition 
is necessary to determine the scattering parameters. 

The paper is organized as follows. Firstly, we will 
characterize the electric polarizability of a dielectric rod. 
Then, we will derive the homogenization rudiments nec- 
essary to calculate the effective permittivity of the struc- 
ture. In section IV, the electromagnetic modes supported 
by the rod medium are characterized. In section V, we 
study the reflection problem in a finite slab, assuming 
that the rods are parallel to the interfaces. The derived 
results are compared with full wave numerical simula- 
tions. Finally, in section VI the conclusions are pre- 
sented. 

In this work we assume that the fields are monochro- 
matic with time variation e'^^"*. 



II. POLARIZABILITY OF A DIELECTRIC ROD 

Let us consider a ENG rod with radius R and (rela- 
tive) permittivity e = e{ui). The rod is oriented along 
the z-direction. In this section, we calculate the com- 
ponent azz of the electric polarizability tensor (actually, 
since the rod is infinite along z, we are going to calculate 
the polarizability per unit of lenght). To this end, we 
consider that a plane wave with magnetic field along the 
y-direction illuminates the dielectric rod. The wave vec- 
tor of the incident wave is k = (fc^, 0, kz), where k.k = 
and (3 = w/c is the free-space wave number. The incident 
electric field along the z-direction is: 

Ei'''' = Eoe-^'''^''e-^''''' (1) 

OO 

n= — oo 

where J„ is the J-Bessel Function of first kind and order 
n, kpfi = — k1, and (r, (f) form a system of cylindri- 
cal coordinates attached to the rod axis. 

The field components along z can be expanded into 



cylindrical harmonics. For example, 

oo 

E anJ\n\ {kp,mr) eJ•"^e-^■'=-^ r<R 
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n= — oo 

(2) 

where a„ and 6„ are the unknown coefficients of the ex- 



pansion. H; 



(2) 



Jn —jYn is the Hankel function of second 



kind and order n, and kp^m = —j^/^z — /3^£. The field 
Hz has a similar expansion. 

The transverse fields E|[ and H|| (projections into the 
xoy plane) can be written in terms of Ez and Hz'. 
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{-jkz'^llEz+jPuz X ViWoHz) (3) 



£ and kp,i = kp^m inside the rod, and 



?7oH|| = -2- {-jPsiUz X V\\Ez-jkz ^wVoHz) (4) 

In the above, % is the free-space impedance, V|| = 
^ Q\ i 

£i = 1 and kp^i = kp^ outside. 

The unknown coefficients can be calculated by impos- 
ing the continuity of the tangential electromagnetic fields 
(i.e. Ez, E^, Hz, H^) at the interface r = R. 

The {z component of the) electric dipole moment (per 
unit of length) is given by. 
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and so the polarizability per unit of length is. 



= (e - 1) 27ri? 



x=y=0 

ao Ji (kp.mR) 



En 



(6) 



This result shows that ao is the unique coefficient re- 
quired to calculate the polarizability. Since Maxwell 
equations are separable in cylindrical coordinates, to cal- 
culate ao it is sufficient to impose the boundary condi- 
tions to the terms associated with the cylindrical har- 
monic n = 0. It can be easily verified that the n = 
term of the magnetic field Hz vanishes (this happens be- 
cause iJ™"^ = and ^ = for the n = harmonic). 
Thus, only the tangential components Ez and H^ have 
non-trivial n = coefficients. Using (2)-(2) in (3)-(4), 
and imposing the continuity of the tangential fields, we 
find that, 

ao Jo {kp,mR) = hoH^o^ [kp^R) + EqJo {kp,oR) (7) 



ao-r^4 ikp,mR) = feoT^i^o^"^ {kp,oR)+-^E^4 (kp^oR) 



kp,o 



kpfi 



(8) 
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where the prime " "' denotes the derivative of a function. 
Solving for ag wc readily obtain: 
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(9) 



Substituting this result in (6), we obtain the desired 
electric polarizability: 
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4 £-1 •■ V ' "'""Jiikp^mR) 
+ ekp^aHl^^ {kp,^R) 



(10) 



In this work, we consider that the radius of the rods 
is always much smaller than the wavelength, or equiva- 
lently, Rkp^ « 1. In these circumstances, (10) simpli- 
fies to, 
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(11) 



where C is the Euler constant. It is important to note 
that because the rods are infinitely long, the polarizabil- 
ity depends not only on the frequency of the incoming 
wave, but also on the wave vector component kz- If e 
approaches —oo the above result reduces to the case of 
perfectly conducting rods studied in [? ]. 

Proceeding similarly, we can obtain the well-known for- 
mula for the electric polarizability (per unit of length) in 
the transverse plane {xoy plane): 
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(12) 



Provided the permittivity of the rods is not too close to 
the resonance e = — 1 and the rods are very thin, the 
polarizability can be neglected in the transverse plane. 



III. EFFECTIVE PERMITTIVITY MODEL 

In the following, we derive an effective permittivity 
model for the medium formed by a periodic array of ENG 
rods. As depicted in Fig. 1, the rods are arranged in a 
square lattice and the spacing between the rods (lattice 
constant) is a. As is well-known, each electromagnetic 
(Floquet) mode in a periodic medium can be associated 
with a wave vector k = {k^, ky,kz). For convenience, we 
define k|| = {kx,ky,Qi). 

To compute the effective permittivity we use the mix- 
ing formula : 



= = 1 = 

£ = I + T ("e -Cint) 

Acell ^ ^ 



(13) 



where Aceii = a?, I is the identity dyadic, the super- 
script "-1" represents the inverse dyadic, and Cint is the 
interaction dyadic calculated in Appendix B. The mix- 
ing formula is derived in Appendix A, and is valid under 
the condition that the dimensions of the cross-section 
are much smaller than the lattice constant, and that 
|k|||a << 2n. As discussed in Appendix A, even though 
(13) reminds Clausius-Mossotti formula, things arc not 
so plain because the lattice has some intrinsic dispersion. 
To keep the readability of the paper the details have been 
moved to Appendix A. 

Next we substitute (11) and (12) in (13). Using (B4), 
we find that the effective permittivity in the transverse 
plane is 



— ^vv — 1 "I" 



J_£+l _ 1 
fv e-1 



(14) 



where fv = nR^/Acdi is the volume fraction of the rods. 
Provided the permittivity of the rods does not satisfy e « 
— 1 and the rods are very thin, we can assume that e^x = 
Cyy ~ 1. For simplicity, we shall assume this situation in 
the rest of the paper. On the other hand, from (13) it is 
clear that, 



= 1 + 



1 



i-cell az 



(15) 
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where Cint,zz is given by (B7). So, after further simplifi- 
cations, we obtain that: 



1 



(£-l)/v 



(16) 



where (3p is the plasma wave number defined consistently 
with the results of [? ] for perfectly conducting wires: 



27r 
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111(2^) +0-5275 



(17) 



Note that in general e is a function of frequency. Formula 
(16) gives the effective permittivity of an array of diluted 
ENG rods. The first important observation is that the 
medium is spatially dispersive. Indeed the permittivity 
depends not only on the frequency fi = oj / c, but also on 
the component of the wave vector parallel to the rods. 
This means that the medium is nonlocal, i.e. in the spa- 
tial domain the electric displacement vector and the elec- 
tric field are related through a spatial convolution rather 
than by a multiplication. Secondly, we note that if the 
permittivity of the rods approaches —00, (16) reduces to 
the formula derived in [? ] for the wire medium, consis- 
tently with the observation made in the introduction of 
this paper. Finally, if we put /? = and assume on-plane 
propagation, i.e. fcz = 0, the effective permittivity sim- 
plifies to Ezz = 1 + (e — 1) /y, which is the exact formula 
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in the static the limit for non-dispersive dielectrics [? ]. 
Thus, very interestingly, formula (16) is in a certain sense 
the average of these two limit situations. Note that even 
though in this work our main interest is the analysis of 
ENG rods, the proposed model is also valid for dielectrics 
with positive real part of the permittivity. In the next 
sections, we characterize the electromagnetic modes sup- 
ported by the rod medium and validate the model with 
numerical simulations. 





IV. CHARACTERIZATION OF THE 
ELECTROMAGNETIC MODES 



FIG. 2: Reflection of a plane wave by a semi-infinite rod 
medium. a)The axes of the rods are parallel to the interface, 
b) The axes of the rods are normal to the interface. 



It is evident that the waves in the homogenized 
medium can be decomposed into transverse electric (TE) 
modes (electric field is normal to the axes of the ENG 
rods), and transverse magnetic (TM) modes (magnetic 
field is normal to the axes of the ENG rods). As ex- 
plained in the previous section, we shall assume in this 
paper that exx = £yy ~ 1 (i-e. that the rods are very thin, 
and that the permittivity of rods is not close to e « — 1). 
Within this approximation, the TE-modes do not inter- 
act with the rods, and thus their dispersion characteristic 



IS, 



(18) 



and the associated average electric field is (the propaga- 
tion factor e^^^ '' is implicit) 



SI 



(19) 



On the other hand, the TM-modes satisfy the character- 
istic equation: 



iP' - ki) 



(20) 



The above equation cannot be solved explicitly as a func- 
tion of (3, because the permittivity of the rods is itself a 
function of [3. The corresponding average electric field is 
(for fc, ^ 0): 



E' 



TM 



SI 



/3 /32e, 



fc2 (3 



(21) 



The associated magnetic field can be calculated using 
(A2). 

To understand better the nature of the TM-modes, 
next we study a reflection problem. Let us consider a 
semi-infinite rod medium illuminated from the air side 
with a plane wave. We analyze two different geometries, 
as depicted in Fig. 2. 

Firstly, let us suppose that the axes of the rods are par- 
allel to the interface x = (Fig. 2a). The incident wave 

vector is k™'^ = {—jlo, ky,kz) with 70 = y/fc^ + — (3'^. 

It is well-known that the component of the incident wave 
vector parallel to the interface, (0, fcjy,/;^), is preserved. 



In this case only one TM-mode is excited in the artificial 
medium (besides the TE-mode). Indeed, from (20) the 
component kx of the wave vector inside the rod medium 
is given by: 



rod 



(22) 



Notice that the right-hand side of the above equation 
only depends on the geometry and parameters of the 
medium, on the wave number (3 of the incident wave, 
and on the components of the incident wave vector that 
are preserved (0, fcy, fc^). 

Next suppose that the rods arc normal to the interface 
z — Q (Fig. 2h). The incident wave vector is now k™"^ = 

{kx,ky, ~j^o) with 7o = ^Jk^~+k^~~P^ . The interesting 

thing is that for this configuration two TM-modes can 
be excited inside the rod medium. Indeed, since k|| = 
{kx,ky,0) to find the excited electromagnetic modes one 



needs to solve (20) for kz 
using (16), show that: 



kl^ (3^ 



Straightforward calculations. 
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^j[(3l + kl-pif +A(3lkl\ (23) 



where we defined the parameter (3c = Pd^) as. 



Pl = 



01 



{e{LS)-\)jy 



(24) 



Note that provided the permittivity of the rods is less 
than the permittivity of the host medium, /3c is a positive 
real number (with the same unities as /3; for simplicity 
the rods are assumed lossless, otherwise /?c becomes a 
complex number). Also, (3c is in general a function of 
frequency since e also is. From (23) it is seen that there 
are two different solutions for kz-, and hence two TM- 
modes, besides the TE-mode, can propagate inside the 
artificial medium. This phenomenon is a manifestation 
of spatial dispersion, and is also characteristic of the wire 
medium [? ] . We also note that the average electric field 
for both TM-modes is calculated using (21). 
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FIG. 3: Plot of kz as a function of normalized frequency /3a, 
for R = 0.05a and fe|| = 0. The permittivity e follows a 
Drude type model (see the text). The solid line represents 
the "exact" result, whereas the dashed line represents the 
results calculated with the proposed model. 



To illustrate the discussion, we plot in Fig. 3, kz as 
a function of normalized frequency /3a for the parame- 
ters R — 0.05a, =0, assuming that the permittivity 
follows the Drude model e = 1 — /3^//3^ with (normal- 
ized) plasma wave number /3ma = 12.0. The dashed line 
curve represents the results calculated using (23). The 
solid hue curve corresponds to the data calculated by 
substituting (10) in (15) (with no approximations) and 
calculating Cn-it,zz using (B5) with r = r' = 0, but with- 
out making any assumptions with respect to k||a being 
negligible. Hence, the permittivity becomes a function 
of not only (3 and kz, but also of the other components 
of the wave vector. The permittivity function obtained 
in this way is substituted in (20), and the corresponding 
equation is solved numerically. A similar procedure was 
used in [? ? ], and so further details are omitted here. 
These results can be regarded as "exact" within the thin 
rod approximation. As seen in Fig. 3, the agreement 
for the first TM-mode is always excellent. On the other 
hand, the second TM-mode is not so accurately predicted 
for relatively small wavelengths. The reason is twofold. 
The first reason is that for larger frequencies the long 
wavelength limit approximation is not so accurate. The 
second reason will be discussed later. Fig. 3 shows that 
for = one of the electromagnetic modes propagates 
with the speed of light. It is also important to refer that 
for very low frequencies one of the TM-modes is cut- 
off (complex imaginary propagation constant). This is 
because the ENG rods behave as perfect conductors in 
the static limit (when the permittivity follows the Drude 
model) . 

To give further insight about the TM-modes, let us 
study different limit situations. First, suppose that at 
some frequency (3c{ijj) « /3p, i.e. the permittivity of 
the rods is very large in absolute value. In this case, 
(23) reveals that one of the modes has the dispersion 
characteristic k1 = 0^ , and that the other mode has the 



dispersion k1 = 0^ — 0^ — k'^y The former mode can be 
readily identified with the well-known transverse electro- 
magnetic (TEM) dispersionless mode of the wire medium 
(perfectly conducting wires), while the latter is the TM- 
mode of the wire medium. 

Consider now the case k|| « 0, i.e. paraxial incidence. 
Using a Taylor expansion we obtain: 
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(25) 



The above formula shows that if either pc{(^) << Pp or 
/3c{w) >> (3p, one of the modes becomes dispersionless 
with respect to k|| (i.e. the coefhcient associated with 
k^ vanishes). The former case was already discussed. 
As to the latter case, the pertinent mode has dispersion 
kl « P'^ + (3^. But this imphes that kz > P and thus 
this mode is a surface wave guided along the rods. In 
[? ] it was proved that a ENG rod is able to support 
tightly bounded surface modes that propagate electro- 
magnetic energy with subwavelength beam radius. For 
metallic materials the surface modes are surface plasmon 
polaritons. It was shown that the energy becomes more 
confined to the vicinity of the dielectric waveguide when 
the effective index of refraction n^ff — kzjP increases. 
This important result justifies why in the rod medium 
one of the TM-modes becomes independent of k||. In 
fact, when /3c (t^) >> Pp each guided mode is confined 
to a small vicinity of the respective ENG rod, there is 
no interaction or coupling between the rods, and conse- 
quently one of the TM-modes becomes dispersionless. As 
referred before, when /3c{uj) << Pp there is also a quasi- 
TEM dispersionless mode. However this mode is qual- 
itatively very different from the mode that arises when 
/3c(w) >> (3p. Indeed, while in the latter case the en- 
ergy is propagated tightly bounded to the ENG rods, in 
the TEM mode case the field energy is distributed more 
or less uniformly by the whole volume. As discussed in 
the introduction, the dispersionless modes may be used 
to canalize the electromagnetic radiation through the rod 
medium and achieve sub-wavelength imaging. A detailed 
analysis of this topic is out of the scope of the present 
paper, and will be reported elsewhere. We also refer that 
the other TM-mode, still assuming that Pc{oj) » Pp, 
has to a first approximation the dispersion characteristic 
kz~P'^ — k'^y i.e. approximately the same dispersion as 
the TE-mode. 

In Figs. 4 and 5, kz is plotted as a function of k^ for 
R = 0.05a and several different values of the normalized 
frequency (3a (the permittivity of the rods follows the 
same Drude model as before). For convenience we show 
the results in two different figures. For very long wave- 
lengths Pd^^) << Pp, and consequently one of the modes 
is cut-off. The dispersion of the other mode (shown in 
Fig. 4) becomes increasingly flat as the frequency (and 
consequently /3c) decreases. Note that from (24) and as- 
suming a Drude type model, (3c increases monotonically 
with the frequency. 
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FIG. 4: Plot of fcz as a function of k^a for R — 0.05a and ky — 
0. Contours /3a = 0.1, 0.3, 0.5, 0.7, 0.9 for the TM-modes. The 
permittivity e follows a Drude type model (see the text). 




FIG. 5: Plot of fcz as a function of kxCi for R = 0.05a and 
ky = 0. Contours /3a = 1.0, 1.2, 1.4, 1.6, 1.8, 2.0 for the TM- 
modes. The permittivity e follows a Drude type model (see 
the text). 



At some point, as the frequency increases, the TM- 
mode that was cut-off starts to propagate. In this case, 
(Fig. 5), for each fixed frequency there are two differ- 
ent contours, i.e. two propagating TM-modes. As the 
frequency increases even more, /3c becomes comparable 
or larger than Pp, and the band structure of one of the 
TM-modes becomes practically flat, consistently with the 
fact that the energy propagated by this mode is tightly 
confined to the vicinity of the ENG rods. 

In Figs. 4 and 5 it can also be seen that around fc|| « 
the wave normal contours of one of the modes are to 
some extent hyperbolic. In fact, using (25) it can be 
easily checked that one of the modes has dispersion char- 
acteristic such that the signs of the coefficients associated 
with and fc^ are symmetric. This property is more im- 
portant for frequencies such that /3c(w) ~ Pp. Note that 
if the medium was anisotropic, with no spatial disper- 
sion, and negative permittivity along the z-direction the 
contours would also be hyperbolic. It is well-known that 



hyperbolic contours may originate negative refraction at 
an interface. 

To conclude this section, we will discuss the scope of 
application of the permittivity model (16). We remember 
that the results were derived for thin rods, R << a, and 
under the assumption that |fcp,o|fi << tt and /3a << tt. In 
general the modes that propagate in the long wavelength 
limit satisfy the previous conditions without problems. 
However, there is one exception at the problem at hand. 
In fact, when /3c(w) >> Pp the radial constant kp Q be- 
comes complex imaginary for the TM-mode associated 
with the surface mode (surface plasmon) . As discussed in 
[? ] , when the beam radius is subwavelength the effective 
index of refraction n^ff — kz/P becomes very large, and 
in that case the condition |A;p,o|a << tt may not be ob- 
served. This situation affects the accuracy of our model 
when Pc{(^) >> Pp. Indeed, the error in the TM-mode 
associated with the surface mode becomes non-negligible 
in this situation. The other TM-mode is still accurately 
predicted. This result justifies deterioration of the agree- 
ment in Fig. 3, as the frequency (and consequently, for 
the Drude Model, also Pc) increases. 

Fortunately, it is easy to solve this problem. In fact, 
when /3c (w) » Pp the dispersion characteristic of the 
pertinent TM-mode is essentially the same as the disper- 
sion characteristic of the guided mode supported by a 
single ENG rod. This dispersion characteristic is deter- 
mined in [? ], and is equivalent to the condition a~} = 0, 
where aj^ is given by (10). Thus, to summarize our find- 
ings, the TM-modes can be accurately calculated using 
(23), except when /3c (cj) » Pp which yields less accurate 
results for the mode with higher . In this case the cor- 
responding TM-mode is dispersionless, and follows the 
same characteristic as the guided mode supported by a 
single dielectric rod. 



V. THE REFLECTION PROBLEM WITH RODS 
PARALLEL TO THE INTERFACES 

To further validate the proposed permittivity model, 
we will study the reflection of electromagnetic waves by 
a rod medium slab with finite thickness. We will suppose 
that the rods are parallel to the interface (see Fig. 6). 
The slab consists of A''^ layers of rods. The structure is 
periodic in y and z, and the dielectric rods stand in free- 
space. The interfaces x = x'j^ and x = x'^ = x'j^ + N^a 
are represented by the dashed lines (since the rods stand 
in air the definition of the interfaces is a bit ambiguous; 
this will be discussed below with more detail). The rods 
in the leftmost layer are in the plane x — xq. 

As discussed in the previous section, when the rods are 
parallel to the interface, equation (22) has only one solu- 
tion for k^. For simplicity, we will restrict our attention 
to the case in which either k^ = (Fig. 6a) or ky — 
(Fig. 6b). For these particular geometries, an incident 
plane wave polarized as depicted in Fig. 6 can only excite 
the TM-mode inside the rod medium. As is well-known 
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[? ], in the general case where ky and kz are simultane- 
ously different from zero, both the TM- and TE-modes 
are excited (the medium is birefringent). Apart from the 
more heavy notation, the general case poses no additional 
difficulties. 
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FIG. 6: A plane wave illuminates a slab of the rod medium 
(Nl ~ 2). The rods are parallel to the plane a; = 0. a) kz = 0. 
b) ky = 0. 

Using (A2), (21), and (22) and matching the tangen- 
tial components of the electric and magnetic fields at the 
interfaces, we find that the reflection coefficient referred 
to the plane x — x'j^ is given by: 



P = 



^ tanh {jmd) (7g - 7^) 
2707m + tanh {jrud) (70 + 7m) 
tanh (7^d) (£^7g - 7^,) 
2£7o7m + tanh (7^^) (£2-^2 ^ ^2^) 



where d = N^a is the thickness of the slab, 7™ — j k. 



(26) 
(27) 

■x.rod-) 



fc2 



h fc2 — /32, and (26) corresponds to Fig. 6a) 
0, and (27) corresponds to Fig. 6b) with 



70 = 
with fcj 
ky = 0. 

Next, in order to demonstrate the accuracy of the the- 
oretical results, the analytical model is tested against full 
wave data computed with the periodic moment method 
(MoM) [? ]. In the first example we consider that 
R = 0.05a, and e = —30.0 at (3a = 1.0 (for simplicity, 
losses are neglected) . A plane wave polarized as depicted 
in Fig. 6a illuminates 5 layers of rods {Nl = 5). The am- 
plitude of the reflection coefficient is depicted in Fig. 7 
as a function of kyU. Note that the angle of incidence is 
such that siny; = ky/ f3. For ky > (3 the incident wave is 
evanescent. The solid line represents the MoM full wave 
data. The dashed line represents the results computed 
using the proposed permittivity model (formula (26)). It 
is seen that that the agreement between the two sets of 
data is good. Similar agreement is obtained for the phase 
of the reflection coefficient and for the transmission co- 
efficient. The previous results also demonstrate that the 
homogenization model is useful to study not only incident 
propagating plane waves, but also part of the evanescent 
spectrum. 

As noted before, since the rods stand in free-space the 
position of the interfaces and thickness of the slab are a 
bit ambiguous. Notice that the thickness of the homog- 
enized slab was taken equal d = N^a, apparently with 
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FIG. 7: Reflection coefficient as a function of the wave vector 
component parallel to the interface. The slab consists of Nl ~ 
5 layers, and the rods have R = 0.05a, and e = — 30 at /3a = 
1.0 The solid line represents the full wave MoM data, and the 
dashed line represents the data obtained using the analytical 
model. 
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FIG. 8: Reflection coefficient as a function of kyu for Nl ~ 1 
layer, R = 0.05a, and e = —30 at /3a — 1.0. The solid and 
dashed lines are defined as in Fig. 7. 



good results. Next, to test if this choice still yields accu- 
rate results for very thin slabs, the reflection coefficient 
is computed for the same structure, except that now the 
slab has only one layer of rods {Nl = 1). The reflec- 
tion coefficient is depicted in Fig. 8 and has a peak at 
kyU = 1.0, which corresponds to the transition between 
propagating waves and evanescent waves. As seen, even 
though the slab is so thin the agreement is still remark- 
ably good. This is a bit a surprising, because for such a 
thin slabs one would expect that the interface effects and 
granularity of the artificial medium would prohibit the 
homogenization of the structure using the bulk medium 
average fields. 

In the next example, we study what happens if the an- 
gle of incidence (p is kept constant {(p = 45[deg]), and the 
frequency is varied. Now R = O.Ola, Nl — 5, and the 
permittivity follows the Drude model e — l~ 0^ with 
plasma wave number (3mO- = 12.0 or = 80.0. The 
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calculated results are shown in Fig. 9. For relatively low 
frequencies the two sets of data agree very well, but as 
the frequency increases the agreement progressively de- 
teriorates, since the long wavelength limit approximation 
is no longer valid. Notice that for relatively low frequen- 
cies the medium blocks the incident radiation because 
the rods effectively behave as perfectly conducting wires. 
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FIG. 9: Reflection coefficient as a function of the free-space 
wave number /3a for R = O.Ola, and Nl = 5 layers. The 
permittivity of the rods follows a Drude-type model (see the 
text). The solid and dashed lines are defined as in Fig. 7. 
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FIG. 10: Reflection coefficient as a function of kztt for Nl = 5 
layer, R = 0.05a, and e = —30 at /3a = 1.0. The solid and 
dashed lines are defined as in Fig. 7. 

So far the wave vector of the incident wave was al- 
ways perpendicular to the axes of the ENG rods, and so 
the effects of spatial dispersion were hidden. In the last 
example, we consider the very different propagation sce- 
nario depicted in Fig. 6b. The parameters of the rods 
are R = 0.05a, e = —30, and Nl = 5. The reflection 
coefficient for /3a = 1.0 is shown in Fig. 10 as a function 
of kza. Note that the angle of incidence 6 for propagat- 
ing waves satisfies sm9 = fcz//3. Fig. 10 shows that the 
agreement between the numerical and analytical results 
is still quite satisfactory, except near the transition be- 
tween the propagating waves and the evanescent waves 



(fc,a« 1.0). 

It is important to refer that the configurations stud- 
ied in this section (rods parallel to the interface) are not 
appropriate to achieve sub- wavelength imaging. For that 
application the rods must be perpendicular to the inter- 
face, as shown in Fig. 2b. As discussed in section IV, 
for this geometry three modes are excited in the (homog- 
enized) rod medium, and thus an additional boundary 
condition is necessary to solve the scattering problem. 
The same situation occurs in the wire medium [? ] . The 
analysis of this problem is out of the scope of this work 
and will be presented elsewhere. 



VI. CONCLUSION 

We studied the electrodynamics of a periodic array of 
thin ENG rods. Using the polarizability of a single rod 
and integral equation methods, we derived new nonlo- 
cal permittivity model for the artificial medium that ac- 
curately describes the propagation of waves in the long 
wavelength limit. We discussed the effects of spatial dis- 
persion in the context of the reflection problem. It was 
proved that when the rods are parallel to interface only 
two modes (TE and TM) can be excited in the artifi- 
cial medium. However, when the axes of the rods are 
normal to the interface, two TM-modes, besides the TE- 
mode, can be excited in the medium, as a manifestation 
of spatial dispersion. It was demonstrated that the wave 
normal contours of one of the TM-modes are intrinsically 
hyperbolic. It was proved that the rod medium supports 
dispersionless modes that propagate along the axes of the 
rods, and it was speculated that this important prop- 
erty may allow sub-wavelength imaging of electromag- 
netic waves at the infrared and optical domains, using 
the idea proposed in [? ] . It was shown that the energy 
of the dispersionless modes can be loosely bounded to 
the ENG rods (in this case the wave is essentially trans- 
verse electromagnetic), or alternatively tightly confined 
to a vicinity of the rods. In the latter case, the modes 
have approximately the same dispersion as the guided 
surface mode supported by an individual rod. The re- 
flection problem was investigated in detail for the case 
where rods are parallel to the interfaces. The developed 
theory was successively tested against full wave data cal- 
culated with the MoM. 
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APPENDIX A: DERIVATION OF THE MIXING 
FORMULA 

In this Appendix we derive the mixing formula (13) 
used to homogenize the rod medium. 

Let us consider a generic electromagnetic Floquet 
mode (E, H), associated with the wave vector k = 
{kx, ky,kz) and the wave number (3 = co/c, i.e. the fields 
satisfy Maxwell-Equations and (E, H) exp (j'k.r) is a pe- 
riodic function in the lattice. The average electric field 
is defined as, 

Eav = y" E (r) '^ '^d^r (Al) 



and Hav is defined similarly. In the above, fl represents 
the unit cell of the periodic medium. Prom Maxwell- 
Equations it can be proved that: 



-k X Eav + /3 %Hav = 



/? ( Eav + — ) + k X r?oHa 







(A2) 
(A3) 



where the generalized polarization vector is given by, 
P 

£0 Vcell 



/ (e- l)Ee+^ '' ''rf^r 

cell J 

n 



(A4) 



Straightforward calculations show that these relations 
imply that: 



Ea 



pH-kk P 

fc2 -/32 



(A5) 



The above equations are exact and completely general. 
Next we will apply these results to the rod medium under 
study. To begin with, we note that since the crystal is 
invariant to translations along the z-direction the fields 
depend on z as exp(— jTc^z). Using standard Green func- 
tion methods [? ] , it can be proved that the electric field 
has the following integral representation, 

E (r) = (£ - 1) J (r| r') .E (r') d^r' (A6) 

s 

where the primed and unprimed coordinates represent 
the source and observation points, respectively, S = 
{ {x', y' ,0) : + < i?2| is the cross-section of the di- 
electric rod in the unit cell, and the Green function dyadic 
is defined by. 



^p=[^+p2 - j''^'^^) i^W - jk,u,)j % (A7) 

In the above, $p = $p(r|r') is the dynamic potential 
created by a phase-shifted array of line sources, 

+ pHp = - ^ J (r|| - rj| - n) e-^''- (■-•') (A8) 



where S is the Dirac function, I = {ii,i2) is a double 
index of integers, rj = a{ii,i2;0) is a lattice point, and 
r|| = {x,y,0). Thus $p is intrinsically two-dimensional, 

depending on z and z' as g^J'''- (-'^-' ). Furthermore, it is 
obvious that the Green function only depends on the rel- 
ative position r|| — rJi • Note that the Green function can 
be written as a superimposition of the potentials created 
by the line sources: 



$p = (r|| -r'li -ri)e 



(A9) 



where $o = e^^''=-(^^^')i7^'^^ {kp^ |r|| -r'|||) /4j is the 
potential created by a line source placed at rj|, i.e. the 
solution of (A8) when the summation in the right-hand 
side is restricted to the index 1 = 0. Physically, (A6) 
establishes that the electric field at some point of space 
is the superimposition of the fields radiated by all the 
dielectric rods of the lattice. 

The Green potential can also be written as a Fourier 
series since it is a (pseudo) periodic function of the wave 
vector: 



$p(r|r') 



1 



-jkj.(r-r') 



Accii ^ kj.kj 



/32 



(AlO) 



where Aceii = a^, J = (ji, j^) is a double index of inte- 
gers, kj = k -I- kjj, and k^ = 27r/a (ji, .72, 0). 

Now that the necessary theoretical formalism was in- 
troduced, we are ready to study the homogenization 
problem in the rod medium. To begin with, we note that 
from (A9) and (AlO) the Green potential is singular in 
the spatial domain, i.e. when r|| — r|| ~ (source region), 

as well as in the spectral domain, i.e. when k.k « 0'^ 

(long wavelength limit). Since the integral (A6) is de- 
fined over the source region and we want to study the 
electromagnetic modes that propagate in the long wave- 
length limit, it is convenient to single out the terms that 
make the Green function singular and decompose it as 
follows: 



$p = $0 + 



-jk.(r-r') 



Acell k.k-/J2 



reg 



(All) 



where $reg, which is defined implicitly by the above for- 
mula, is a regular function both in the spatial domain 
(source region r|| — rj| w 0) and in the spectral do- 
main (long wavelength limit). Using this decomposition 
in (A6) we find that: 



E(r) = 



Eave-^'^-'- + 



(£-l)y /?2Go(r|r').ErfV 

s 



+ 



(e-l)y f3^Greg{r\r').Ed^r' 
s 



(A12) 



where Eav is the average field in the crystal, and Gg and 
Greg are defined as Gp, except that the Green potential 
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is replaced by <I>o and ^reg, respectively. To obtain the 
above result we used (A5). Next, we use the fact that the 
ENG rods are assumed to be very thin {R/a << 1), and 
that we want to investigate the electrodynamics of modes 
that propagate in the long wavelength limit (|k|| |a << 27r 

and (3a « 2w). Since the dyadic Greg is regular in both 

the spatial and spectral domains, it is Icgitim to write 
(putting r|| = rJi = and k|| = {kx,ky,0) = 0), 

%eg (r| r'; k, /?) ^ %eg iz\z'; k,,P) (A13) 

being the formula valid in the source region. For conve- 
nience, we introduce the following interaction dyadic: 

^i^t=P^%egiO;k„p) (A14) 
Then, it is clear from (A12) that, 

E(r)« (^Eav+a„t.^)e-^'=^^ 

+ (e - 1) y" p^^o (r| r') .E dh' (A15) 
s 

provided the observation point r is near the source region 

and |k|||a << 27r. In above, we introduced the electric 
dipole moment (per unit of length), p, of the dielectric 
rod in the unit cell. Now, the key result is that (A15) 
is formally equivalent to the (integral) equation obtained 
when a single rod is illuminated by a plane wave with 

electric field with amplitude E™'^ = ^Eav + Cint-^) and 

the same wave vector component k^ along the z-direction. 
In other words, when the rod in the unit cell stands alone 
in free-space and is illuminated with the above defined 
plane wave, the total electric field also satisfies (to a first 
approximation) (A15) in the source region. But this re- 
markable result implies that: 

^ =Se. fEav+Oinf-) (A16) 

where cie is the electric polarizability tensor for a single 
rod. The term inside brackets in the right-hand side can 
be regarded as the local field that polarizes a single rod 
embedded in the dielectric crystal. Within the thin rod 
condition and for long wavelengths, the above solution is 
exact. 

We are now ready to calculate the effective permit- 
tivity dyadic. Since the (macroscopic) electric displace- 
ment vector D is given by D = eoEav + p/Aceii, and 
the effective permittivity tensor must guarantee that 
D = eoe.Eav, we conclude that the effective permittivity 
of the rod medium is given by the mixing formula (13). 
Note that (13) is completely general and is valid indepen- 
dently of the specific geometry of the transverse section 
of the rod. 

At this point it is appropriate to compare (13) with 
the classic homogenization approach. It is striking that 



(13) reminds Clausius-Mossotti formula [? ? ]. In- 
deed, if we could identify the interaction dyadic Cint with 
l/2Aceii the formulas would be the same (note also that 
the rods arc arranged in a square lattice). In Appendix 
B we calculate the interaction dyadic in closed-form us- 
ing the static limit approximation. Equation (B4) shows 
that the interaction dyadic is different from 1 /2Accii only 
along the z-direction. This important result shows that 
Clausius-Mossotti formula is not valid for media with 
cylindrical inclusions. More specifically it fails to predict 
the effective permittivity along the direction in which the 
crystal is uniform. This is an indirect manifestation of 
spatial dispersion. 

We also mention that the interaction dyadic defined 
by (A14) is not equivalent to the dynamic interaction 
constant defined in other works (see for example [? ]). 
Indeed, in our definition we extracted the singularities 
in both the spatial and spectral domains, while other 
works usually only extract the singularity in the spatial 
domain. It is clear from our previous analysis that it 
is the singularity in the spectral domain that indirectly 
defines the relation between the local field that polarizes 
the rod and the average field. 



APPENDIX B: CALCULATION OF THE 
INTERACTION DYADIC 

Here we calculate the interaction dyadic defined by 
(A14). It can be written as: 

Cint = (V|| -jfc^U.) (V|| -jfc.U,)) ^reg (Bl) 

where the right-hand side of the expression is evaluated 
at r = r' = and k|| = 0. Prom (A8) and (All) it is 
clear that: 

V'^^reg + P'^^reg = 

Putting r = r' = and k = in the above equation, and 
letting p approach zero (static limit), we find that: 

V^<^reg (r = r' = 0; k = 0)1 = (^3) 

-f^cell 

Moreover, because of the symmetry of the square lattice 
it is evident that if r = r' = and k = the following 

relations hold, qJ^q^^. — if z ^ j, ^ q^s"" = 0, and 

^ Q^i" = ^ Qyi" ■ So using (B3), we conclude that in the 
static limit (k = and /3 = 0) the interaction dyadic is 
given by: 
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The above result is consistent in the xoy plane with the 
(two dimensional version of the) Clausius-Mossotti for- 
mula. However, along the 2;-direction, maybe a bit sur- 
prisingly, the interaction constant vanishes in the static 
limit. Next, we will estimate Cint,zz in the dynamic case. 
Prom (Bl), we have that: 



Using (A9) and (All), and putting r = r' 
0, we obtain that: 

a„t,z. = - kl) J2 ^^^o'^ (fcp.o |ri|) 



(B5) 
and kii = 



1 



»-cell 



(B6) 



The series in the right-hand side was evaluated in [? ]. 
Using the results of [? ], and assuming that kp^a « n, 
we obtain: 



int, zz 





h-ln( 






2w \ 


V 4^ / 



+ 



1 V A 



12 ' ^ Trlnl e^'rl"! - 1 

n=l ' ' / 



where C is the Euler constant. 



C 



(B7) 



